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When  external  stresses  in  a  system  -  physical,  social  or  virtual  -  are  relieved  through  im¬ 
pulsive  events,  it  is  natural  to  focus  on  the  attributes  of  separate  avalanches.  However, 
during  the  quiescent  periods  in  between,  stresses  may  be  relieved  through  competing  pro¬ 
cesses,  such  as  slowly  flowing  water  between  earthquakes  or  thermally  activated  disloca¬ 
tion  flow  between  plastic  bursts.  Such  unassuming,  smooth  responses  can  have  dramatic 
effects  on  the  avalanche  properties.  Our  thorough  experimental  investigation  of  slowly 
compressed  Ni  microcrystals,  covering  three  orders  of  magnitude  in  nominal  strain-rate, 
exhibits  unconventional  quasi -periodic  avalanche  bursts  and  higher  critical  exponents  as 
the  strain  rate  is  decreased.  Our  analytic  and  computational  study,  naturally  extend¬ 
ing  dislocation  avalanche  modeling  to  incorporate  competing  processes,  reveals  the  emer- 
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gence  of  avalanche  oscillator  scaling  behavior,  a  novel  critical  state  manifesting  through 
self-organized  oscillatory  approaches  toward  a  critical  depinning  point.  We  demonstrate 
that  the  predictions  of  our  theory  are  consistently  manifested  in  our  experiments. 

Physical  systems  under  slowly  increasing  stress  may  respond  through  abrupt  events.  Such 
jumps  in  observable  quantities  are  abundant  (7),  from  Barkhausen  noise  (2)  to  plastic  flow  (3-5) 
and  earthquakes  (6).  Even  though  events  appear  randomly  sized  and  placed,  the  statistical 
properties  of  avalanches  are  universal,  defining  well  understood  non-equilibrium  universality 
classes:  The  main  unifying  concept  is  the  depinning  of  an  interface  under  an  external  field  (7,7). 
An  implicit  assumption  underlying  these  concepts  is  that  all  other  coexisting  physical  processes 
are  either  too  fast  and  thus  average  out,  or  too  slow  rendering  a  static  approximation  valid. 
However,  the  latter  assumption  is  not  always  true  if  the  slow  processes  rearrange  the  pinning 
landscape  at  rates  comparable  to  the  external  field  driving  rates.  For  as  the  fast  avalanches  are 
scale  invariant,  the  whole  timeseries,  including  the  waiting  intervals  between  the  fast  events,  is 
also  scale  invariant  (8).  It  is  there  within  the  waiting  intervals  that  a  slow  restructuring  of  the 
pinning  field  can  thrive  and  alter  universal  predictions  as  observed  by  Jagla  (9). 

While  intermittent  plastic  flow  is  well  known  (10),  only  recently  was  it  shown  as  statistically 
akin  to  universal  mean-field  avalanche  behavior  in  the  quasistatic  limit.  Investigations  of  the 
phenomenon  utilized  a  wide  variety  of  techniques,  including  acoustic  emission  from  deforming 
ice  (3),  high  resolution  extensometry  of  tensile  strained  Cu  (11)  and  microcrystal  compression 
tests  for  FCC  and  BCC  crystals,  including  Ni  (4),  A1  (72),  Mo  (13)  and  FiF  (14).  However, 
most  of  these  studies  covered  only  a  narrow  range  of  nominal  high  strain  rates.  Preliminary  evi¬ 
dence  that  suggests  a  more  complex  physical  picture,  was  discussed  by  some  of  us  in  Ref.  (75), 
where  a  rate  dependence  of  the  cumulative  strain  event  size  distributions  was  observed.  In  our 
experiments,  Ni  microcrystals  of  comparatively  large  dimensions,  having  diameters  between 
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18  and  30  //rn,  were  uniaxially  compressed  (15).  By  controlling  the  applied  external  stress  to 
maintain  a  nominal  strain  rate  and  by  detecting  slip  with  extremely  sensitive  extensometry,  we 
track  the  crystal  displacements  in  time.  In  order  to  study  the  rate  dependence,  we  perform  our 
experiments  at  three  different  nominal  strain  rates  (10_4/s,  10_5/s,  10_6/s).  For  each  sample,  the 
timeseries  of  the  displacement  time  derivative  is  filtered  using  optimal  Wiener  filtering  methods 
adapted  for  avalanche  timeseries  (16),  and  avalanche  events  are  appropriately  defined  without 
using  thresholding. 

As  plastic  deformation  proceeds  in  the  micropillars,  the  dislocation  ensemble  evolves  at  dif¬ 
ferent  time  scales.  The  most  apparent  activity  is  associated  with  fast  glide  processes  which 
produce  stochastic  bursts  of  deformation.  Concurrently  and  in  between  these  events,  other  less 
observable  processes  (cf.  Fig.  1)  contribute  to  collective  slower  relaxations .  Like  glide,  these 
too  are  thermally  activated  processes  accessible  at  these  high  levels  of  stress,  but  having  dif¬ 
ferent  activation  barriers,  for  example:  the  viscoelastic  response  of  the  dislocation  forests  after 
fast  avalanche  strain  bursts,  the  localized  dislocation  climb  motion  in  directions  other  than  the 
glide  plane  under  high  local  stresses,  and  also  the  cross-slip  processes  of  dislocations  shifting 
between  glide  planes  are  slow  processes  compared  to  the  avalanche  events  (17).  They  all  com¬ 
pete  to  minimize  the  far  field  stress  while  changing  the  local  stress  landscape  and  bypassing  the 
glide  process.  They  affect  dislocation  slip,  but  at  a  slower  rate  than  avalanche  glide  (18).  In 
our  experiments,  we  classify  as  “slow  relaxation”  all  the  deformation  that  does  not  belong  to 
avalanches  of  the  scaling  regime.  Using  this  definition,  the  slow  relaxation  fraction  increases 
drastically  at  the  two  slowest  strain  rates.  Thus,  the  rate  dependence  of  the  avalanche  size  dis¬ 
tribution  (Fig.  2(c))  occurs  when  the  nominal  strain-rate  becomes  comparable  to  the  rate  of  the 
slow  relaxation  processes  (cf.  Fig.  1).  Although  the  exact  mechanisms  are  unknown,  one  lo¬ 
calized  reorganization  mechanism  possible  at  these  large  local  stresses  and  low  temperatures 
(0.17Tm  (~  300K))  (18)  could  be  tied  to  newly  discovered  unconventionally  large  cross-slip 
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rates,  as  calculated  for  similar  conditions  to  our  experiments  (19).  For  our  purposes,  and  re¬ 
gardless  of  the  type  of  relaxation  mechanisms,  we  focus  on  the  experimental  fact  that  relaxation 
and  driving  rates  become  comparable.  We  phenomenologically  incorporate  the  slow  relaxation 
into  our  model  in  a  simple  and  intuitive  manner  and  then  a  posteriori  show  that  our  results  are 
independent  of  the  particular  form  of  relaxational  dynamics  (see  SOM). 

The  time  series  of  slip  event  sizes  S,  labelled  by  their  beginning  time,  display  a  striking 
dependence  on  the  driving  rate.  After  we  smooth  the  timeseries  over  a  fixed  window  of  400s 
and  then  rescale  the  time  axis  to  display  comparable  strain  evolution,  a  very  clear  (cf.  Fig.  2(a)) 
quasiperiodic  behavior  emerges  at  the  10~6/s  nominal  rate.  The  emergent  time  period  is  ob¬ 
served  to  display  an  exponential  dependence  on  the  nominal  strain  rate,  while  its  magnitude 
reaches  ~  8  hrs  (for  10~6/s),  consistently  much  larger  than  the  length  we  chose  for  the  fixed 
window  averaging  (cf.  Fig.  1(c)).  This  qualitative  change  in  the  behavior  is  also  reflected  within 
statistical  distributions  of  S:  The  size  distribution  shows  a  power  law  behavior  ( P(S )  ~  S~T) 
for  all  studied  strain  rates  (10”4/s,  10~5/s,  10~6/s),  but  the  value  of  the  power  law  exponent 
drifts  from  ~  1.5  (consistent  with  Refs.  (4,20))  to  a  higher  and  unexpected  value  of  ~  2.0  (cf. 
Fig.  2(c)).  An  analogous  behavior  is  observed  for  the  durations  T  of  the  avalanches  and  their 
correlation  with  the  sizes  (cf.  SOM);  although  there  are  concerns  about  the  accuracy  of  these 
durations,  given  the  response  time  of  the  apparatus  and  the  frequency  of  the  data  recorded  (15). 

Our  explanation  of  the  experimental  data  builds  on  the  model  framework  of  dislocations 
moving  through  a  disordered  landscape  of  forest  dislocations,  on  a  single  slip  plane  under  the 
application  of  shear  stress.  This  is  a  successful  picture  for  avalanches  during  stage  I  plastic¬ 
ity  (20-23)  that  strongly  relies  on  well-understood  models  of  2+1  dimensional  interface  depin- 
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ning  (7).  Our  model  represents  a  minimal  generalization  via  an  added  relaxation  term,  D, 


(i)  (ii) 


where  //  is  the  shear  modulus  of  the  system,  and  e«l.  Here,  0  denotes  the  basic  slip  variable 
of  the  system,  which  is  the  yx-component  (Burgers  vector  along  x)  of  the  plastic  distortion 
tensor,  considering  only  infinite  dislocations  along  z  on  xz  slip  planes  (23).  Part  (I)  of  Eq.  1 
denotes  the  coarse-grained  slow  relaxation  of  edge  dislocations,  assuming  a  rate  D  at  a  fixed 
temperature.  Only  positive  slip  motion  is  considered  to  simplify  our  simulations  (24).  With 
D  ~  exp(— A/T),  we  define  an  effective  rate  of  thermally  activated  processes  (A  is  an  effective 
energy  barrier  to  such  processes)  that  lead  to  slow  relaxation.  For  simplicity,  we  set  the  exponent 
n  —  1,  since  our  conclusions  do  not  qualitatively  depend  upon  it.  The  applied  stress  is  the  xy- 
component  of  the  stress  tensor, 

a(r)  =  (Text  +  <bnt(r)  +  (xhard(r)  =  Met  +  J  d2r'K( r  -  r')0(r')  -  k<j>{ r)  .  (2) 

We  consider  a  stress-controlled  test  in  a  stationary  plastic  regime  (crext  =  Met)  (25),  where  M 
is  a  machine  stiffness,  so  that  c  has  strain-rate  units.  In  this  framework,  the  relative  timescales  of 
the  relaxation  and  stress  rate  increase  are  controlled  by  the  dimensionless  parameter  R  =  D/c. 
Part  (II)  of  Eq.  1  denotes  the  fast  glide  process  which  controls  the  avalanche  dynamics.  Hard¬ 
ening  is  phenomenologically  represented  via  a  coefficient  k  that  controls  the  distance  of  the 
system  from  its  non-equilibrium  depinning  critical  point.  For  clarity,  we  separate  the  relevant 
timescales  by  considering  e«l,  leading  to  infinitely  fast  avalanches  compared  to  the  slow  re¬ 
laxation  process.  <7/  denotes  the  uncorrelated  local  pinning  potential,  due  to  dislocation  forests. 
Finally,  crint  is  chosen  to  contain  the  appropriate  interaction  kernel  K  for  single  slip  straight 
edge  dislocations  (23).  However,  our  main  qualitative  conclusions  are  independent  of  the  ker¬ 
nel,  and  thus  are  equally  applicable  to  other  models  of  avalanches  in  plasticity  (26). 
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The  model  of  Eqs.  1  and  2  is  solved  by  explicit  integration  with  an  Euler  step  until  the 
condition  cr(r)  >  o/(r)  is  satisfied  for  some  point  r,  at  which  point  an  avalanche  is  triggered. 
For  no  relaxation  (D  =  0),  the  avalanches  display  statistics  consistent  with  the  predictions  of  the 
mean-field  theory  of  interface  depinning  (1,5).  As  the  relaxation  I)  increases,  both  the  critical 
exponent  r  for  strain  jump  sizes  S  ( P(S )  ~  S~  T)  and  the  critical  exponent  a  for  event  durations 
T  ( P(T )  ~  T~a),  increase  substantially  (cf.  Fig.  2(b),  SOM).  However,  the  dependence  of 
( S )  on  the  duration  T  (( S )  ~  T1/'7*'2)  remains  unaltered  (1/auz  ~  2),  signifying  that  the 
universality  class  of  the  underlying  non-equilibrium  critical  point  remains  unaltered  while  some 
critical  exponents  increase.  In  the  context  of  mean-field  theory,  remotely  similar  behavior  takes 
place  when  the  driving  rate  c  is  increased  (16)  leading  to  avalanche  overlap,  unrelated  though 
to  our  results  which  hold  when  c  — ■>  0  keeping  R  fixed  (and  >  1). 

The  increase  of  the  exponents  is  accompanied  by  a  quasiperiodic  behavior,  signified  by 
large  slip  events  (cf.  Fig.  2(b)):  If  one  considers  the  average  avalanche  size  in  a  window  (cf. 
Fig.  2(b)),  similar  to  the  experimental  study  but  without  strain  from  relaxation  included,  it 
is  clear  that  the  D  =  0  flat-in-time  profile  is  replaced  by  strongly  oscillating  profiles  in  the 
presence  of  slow  relaxation  (D  >  0).  The  average  avalanche  size  (D  =  0)  is  inversely  related 
to  the  hardening  coefficient  k,  k  ~  (S)~ak  with  crfe  ~  1 .  Thus,  there  is  a  distribution  of 
hardening  coefficients  being  effectively  sampled,  reflecting  local  heterogeneity.  We  assume 
that  such  a  distribution  g(k')  biases  the  integration,  over  all  possible  hardenings  k' ,  of  the  size 
probability  distribution  of  the  D  =  0  model,  leading  to  the  observed  dynamically  integrated 
size  distribution.  That  is,  a  curve  in  Fig.  2(d)  may  be  obtained  as 

roo 

Pint(S)=  /  g(k')P(S,  k')dk'  ,  (3) 

Jo 

For  example,  assuming  that  g(k')  =  const,  for  any  k’  <  hy  and  zero  above  a  limiting  A'0 ,  then 
we  simply  have  at  D/c  1,  I)llt  (S)  =  S~T~CTkV(SklJ(Jk) ,  yielding  a  higher  effective  sizes- 
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exponent  t  =  r  +  for  slow  strain  rates,  with  f  =  5/2  at  mean-field.  However,  typically, 
g(k')  might  have  a  more  complicated  functional  form,  allowing  for  a  large  range  of  possible  fs. 


The  profound  effects  of  slow  rate  processes  within  our  dislocation  model  and  the  comparison 
with  experiments  forces  us  to  ask:  Are  our  findings  general?  To  make  analytical  progress 
toward  an  answer  consider  the  “susceptibility  to  jumping”  p,  defined  as  the  number  of  degrees 
of  freedom  that  slip  after  the  trigger  slip.  On  average  in  our  model  (p)  is  proportional  to  the 
hardening  coefficient  ((p)  oc  1  —  k).  In  traditional  mean-field  models  of  interface  depinning  ( 1 ), 
this  quantity  describes  the  “distance”  of  the  system  from  the  critical  point  and  is  considered 
saturated  to  a  fixed  point  value  after  short-time  transients  (steady  state).  When  p  < Cl,  the 
system  is  far  from  the  critical  point,  while  the  system  is  near  critical  when  p  ~  1.  Numerical 
solutions  to  Eq.  1  verify  that  the  additional  relaxation  process  affects  p  in  an  unusual  way: 
When  an  avalanche  with  size  St  takes  place,  p  instantaneously  decreases  in  a  way  proportional 
to  the  magnitude  of  St,  while  it  increases  linearly  when  there  are  no  avalanche  events.  In  a 
minimal  sense,  we  suggest  that  the  basic  physical  mechanism  behind  the  behavior  of  Eq.  l(with 
c  — >  0  but  R  fixed)  is  given  by  the  behavior  of  the  susceptibility  to  jumping  p,  whose  basic 
characteristics  can  be  described  by  the  following  Markov  process, 


Pt+ i  ~  Pt  =  Apt  =  cd(  1 


(4) 


where  St  is  picked  from  the  mean-field  size  distribution  P(St)  =  NS^2e"  St/Sl>  (where  N  is  a 
normalization  factor  and  S0  =  a/ (1  —  pt))  (1)  and  the  step  cd  shall  be  thought  as  being  propor¬ 
tional  to  R.  The  well-known  avalanche  mean-field  behavior  is  described  by  the  trivial  cd  0 
fixed  point  (analogous  to  higher  experimental  strain  rates).  The  size  of  the  avalanche  at  time 
t,  S't ,  is  a  stochastic  variable  which  mimics  the  avalanche  dynamics  described  by  Eq.  1  (which 
is  in  the  mean-field  universality  class).  When  cd  <C  1,  p  increases  in  small  steps  towards  the 
fixed  point  given  by  p0  =  1  —  ffi  with  average  size  S  =  f  v  vr/(l  —  p0)  (a  being  the  minimum 
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accessible  size)  saturating  at  a  value  similar  to  that  for  the  traditional  depinning  dynamics  (cf. 
Fig.  3  left).  However,  there  is  a  finite  probability  of  a  large  (but  not  infinite!)  avalanche  which 
takes  the  system  far  from  the  fixed  point,  with  A pt  large  and  negative.  If  SS  =  St  —  S,  then 
A pt  =  —cd5S/S  ~  — 1  indicates  the  emergence  of  a  novel  quasiperiodic  behavior  (cf.  Fig.  1(b)) 
showing  large  negative  p— jumps  with  St  being  large  rare  avalanche  events,  much  larger  than 
S.  Now,  p0  is  not  able  to  have  the  essential  defining  role  as  a  steady-state  susceptibility  as  for 
the  traditional  critical  behavior,  and  p  performs  a  Sisyphean  task  constantly  ascending  towards 
p0 .  In  this  way,  the  distribution  of  p  effectively  flattens  as  cd  increases  (low  experimental  strain 
rates)  (cf.  Fig.  3  right),  leading  to  integrated  exponents  (cf.  Eq.  3).  In  a  consistent  manner,  the 
analogous  distribution  for  the  model  in  Eq.  1  flattens  as  R  increases  (cf.  Fig.  3  right  (inset)). 

The  rare  5S  events  scale  with  S0  ~  4,S'2/ (7ra)  and  qualitatively,  the  transition  between  the 
two  regimes  takes  place  when  cd  ~  a7r/(4S)  ~  l/(S/a).  Numerically,  we  get  a  consistent 
“phase  diagram”  line  (cf.  Fig. 3  lower)  which  separates  the  two  qualitative  behaviors.  We  shall 
mention  that  our  model  shares  a  strong  resemblance  to  the  case  of  relaxation  limit  cycle  oscil¬ 
lations  near  a  singular  Hopf  bifurcation  with  stochastic  perturbations  (27),  a  dynamical  system 
similar  to  Sozu,  the  traditional  Japanese  gardening  device. 

Contemporary  observations  have  revealed  novel  and  creative  mechanical  behaviors  of  crys¬ 
tals  in  the  microscale.  Together  with  the  size  effects  (28, 29)  and  the  emergence  of  avalanche 
slip  events  (3,4),  the  importance  of  often  neglected  slow  processes  on  intermittency  has  now 
come  to  light.  Our  study  provides  basic  elements  for  constructing  a  much  desired  accurate  mul¬ 
tiscale  theory  of  microscale  plasticity  (30).  By  incorporating  the  complex  microscopic  slow 
dynamics  of  dislocations  operating  between  the  fast  glide  processes,  it  is  possible  to  capture  the 
complete  interplay  between  local  stress  concentrations  that  lead  to  violent  bursts  and  the  collec¬ 
tive  relaxational  dynamics  of  the  dislocation  environment  over  longer  times.  These  experiments 
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in  the  world  of  “small”  and  the  associated  theory  that  we  developed,  force  us  to  reconsider  our 
understanding  of  the  world  of  “large”,  earthquake  faults  ( 6 ),  sheared  amorphous  solids  (31) 
and  other  depinning  processes,  setting  additional  constraints.  Whenever  the  mechanism  that 
drives  avalanche  dynamics  competes  with  other  slow  coexisting  processes  to  lower  stress,  when 
the  driving  is  slow  enough,  the  nature  of  the  avalanche  dynamics  is  modified  and  gives  rise  to 
critical  quasiperiodic  bursts. 
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Figure  1  -  Complexity  of  dislocation  motion  and  slow  relaxation  during  avalanche  wait¬ 
ing  intervals.  Upper:  A  schematic  toy  demonstration  of  the  complex  processes  that 
dislocation  experience.  Lighter  to  darker  color  scale  indicates  time  evolution.  Under  the 
application  of  stress,  a  dislocation  loop  nucleates  and  attempts  to  grow  but  it  gets  pinned 
on  its  slip  plane(lighter  color).  Depinning  to  a  larger  loop  configuration  corresponds  to  a 
fast  glide-slip  burst,  the  most  common  dislocation  slip.  However,  as  the  loop  grows  a  screw 
dislocation  segment  undergoes  double  cross- slip  to  a  parallel  slip  plane,  in  order  to  bypass 
the  glide  barriers.  On  the  new  plane,  the  dislocation  loop  glides  and,  ultimately  it  under¬ 
goes  dislocation  climb.  Lower  Left:  A  part  of  the  timeseries  of  the  plastic  deformation 
of  a  Ni  sample  at  1CU6  rate.  The  avalanche  phenomenon,  when  close  to  criticality,  does 
not  manifest  itself  only  through  the  rather  fast  and  violent  scale-invariant  bursts  (16),  but 
also  through  unavoidable  long  waiting  times  (8,32)  where  glide  processes  are  inactive;  It 
is  due  to  this  scale-invariant  distribution  of  waiting  times  that  slow,  non-equilibrium  stress 
minimizing  processes  can  thrive,  altering  in  a  fundamental  manner  the  properties  of  the 
critical  avalanche  behavior.  Lower  Right:  The  estimated  percentage  of  strain  accumulated 
in  slow  relaxation  flow  which  does  not  amount  to  large,  observable  strain  bursts  is  shown 
on  the  left  (black):  The  percentage  (relaxation  strain/total  final  strain)  strongly  increases  as 
the  rate  decreases.  Experimental  noise  certainly  contributes  to  the  relaxation  strain  mea¬ 
sured.  On  the  right,  a  non-trivial  quasiperiod  (average  period  of  intermittent  oscillations) 
of  organised  avalanche  behavior  emerges  and  increases  dramatically  as  the  nominal  rate 
decreases. 

Figure  2  -  Experiments  vs.  Theory,  (a)  Average  avalanche  size  as  a  function  of  time  for  different 
strain  rates.  Avalanche  sizes  are  averaged  based  on  their  start  times  in  400s  windows.  Note 
that  the  time  axes  are  rescaled  by  the  nominal  strain  rate,  aligning  the  “strain  scales”.  For 
different,  independently  strained  samples,  clear-  quasi- periodic  avalanche  behavior  emerges 
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as  the  nominal  strain  rate  decreases.  The  period,  as  soon  as  it  becomes  apparent  (0.4nm/s 
sample),  stays  of  similar  size  in  “strain  scale”  while  it  increases  drastically  in  real  time  - 
a  key  prediction  of  the  theory  we  present  in  this  paper,  (b)  The  stick-slip  oscillations  that 
are  observed  experimentally  (a)  are  a  typical  characteristic  of  the  explored  model  of  Eq.  1 . 
Here,  we  fixed  the  value  of  the  relaxation  rate  and  varied  the  actual  strain-rate  (by  modifying 
c),  similarly  to  the  situation  in  the  experiments.  We  should  mention  that  there  is  an  overall 
scale  that  is  freely  modified  in  order  to  match  the  experiments.  This  gives  us  a  unit  of  strain 
being  2  x  1CU6  and  a  fast  timescale  unit  of  0.5s.  We  shall  notice  that  there  is  a  difference  of 
the  simulations’  timeseries  from  the  experimental  in  (a)  in  that  the  detected  events  are  the 
actual  avalanche  events  without  the  distortion  that  would  appeal-  due  to  the  strain  coming 
from  the  slow  relaxation.  Ideally,  this  plot  should  experimentally  emerge  if  the  creep  strain 
is  removed.  This  difference  is  behind  the  overall  scale  difference  of  the  two  plots,  (c) 
Rate  dependence  of  the  critical  exponent  of  the  probability  distribution  for  the  sizes  of  the 
displacement  jumps.  Even  though  the  hardening  coefficient  for  the  different  samples  used 
is  roughly  similar  (±10%),  the  slope  of  the  distribution  increases  as  the  nominal  strain 
rate  decreases.  We  shall  notice  that  as  the  rate  decreases  a  large  collection  of  small,  almost 
invisible  events  appeal-  with  similar  probability  in  the  distribution  (not  shown),  (d)  Similarly 
as  in  (c),  the  decrease  of  c  leads  to  a  ’’rate  dependence”  of  the  size  distribution  exponent  (and 
other  related  exponents),  which  we  explain  in  the  text.  The  exponent  drifts  from  the  known 
~  1.5,  reported  also  in  the  original  Ref.  (23),  to  approximately  ~  2.5,  the  limiting  value 
for  a  ’’flat”  distribution  of  p’s  (cf.  Fig.  3).  As  we  show  in  the  Supplementary  Information, 
it  is  clear  that  the  high  S  bump  near  the  cutoff  fades  away  and  flattens  as  k  increases,  even 
though  the  exponent  stays  virtually  unaltered. 

Figure  3  -  Effective  modeling  of  the  jumping  susceptibility.  Upper  Left:  A  plot  portraying 
the  behavior  of  the  reduced,  analytical  model  Eq.  4.  As  the  rate  q  is  increased,  large  noise 
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causes  p(t)  to  oscillate  between  p  ~  1  and  small  p,  causing  larger  exponents,  according 
to  the  discussion  in  the  main  text.  Upper  Right:  The  probability  distribution  of  p  as  cal¬ 
culated  from  the  solutions  of  the  effective  model  Eq.  4.  The  distribution  becomes  flat  as 
Cfj  increases,  with  p  being  “sampled”  with  equal  probability  during  an  oscillation.  In  the 
inset,  the  probability  distribution  of  (5)50  is  shown  for  the  simulations  of  Eq.  1 ,  calculated 
by  considering  the  timeseries  of  the  events,  regarded  as  instantaneous,  and  then  averaged 
with  a  running  window  of  size  50.  The  simulations  have  the  same  parameters  as  in  (a). 
The  histogram,  shown  in  the  appropriate  scale  (1  —  1  / (S)2  =  p)  becomes  flat,  while  be¬ 
ing  sharply  peaked  for  small  R.  Center:  The  phase  diagram  which  distinguishes  the  novel 
regime  where  the  fluctuations  of  p  are  large  enough  to  cause  quasiperiodic  oscillations  from 
the  traditional  regime  where  p  fluctuates  weakly  around  the  fixed-point  value  po-  The  line 
signifies  the  instability  to  quasiperiodic  behavior  ~  1/S  as  described  in  the  text.  The 
points  label  the  locations  where  simulations  of  Eq.  4  were  performed  and  then,  the  critical 
exponent  t  was  calculated  using  Eq.  3  at  the  points,  with  a  final  Gaussian  interpolation  for 
the  color  background. 
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Figure  1  -  Complexity  of  dislocation  motion  and  creep  during  avalanche  waiting 
intervals 
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Figure  2  -  Experiments  vs.  Theory 
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Figure  3  -  Effective  modeling  of  the  jumping  susceptibility 
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This  PDF  file  includes: 

(A)  Materials  and  Methods 

(B)  SOM  Text 

(C)  Figs.  SI  to  S12 

(A)  Materials  and  Methods 

Experimental:  The  experimental  measurements  were  performed  using  the  methodology  de¬ 
scribed  in  (10).  The  data  are  taken  at  time  resolutions  5,  50  and  500  Hz  for  different  samples, 
depending  on  the  case.  The  nominal  strain  rates  were  10-4,-5,-6/s  with  corresponding  aver¬ 
age  platen  velocities  of  4,  0.4  and  0.04  nm/s,  given  the  dimensions  of  the  pillars.  The  experi- 
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mental  timeseries  were  filtered  using  optimal  Wiener  filtering  methods  optimized  for  studying 
avalanches,  as  described  in  (3).  Optimal  Wiener  filtering  corresponds  to  a  low-pass  filter  that 
has  significant  effects  only  at  short  timescales  which  are  plagued  by  apparatus  problems.  In  a 
similar  fashion,  as  in  Ref.  (3),  we  performed  adequate  tests  in  order  to  confirm  that  the  power- 
laws  and  the  long-time  quasiperiodic  behavior  are  not  related  to  the  filtering  procedure.  The 
experimental  data  is  available  upon  request. 

Theoretical:  In  the  simulations  of  Eq.  1,  during  diffusion,  Euler  time  stepping  is  used  to 
evolve  the  differential  equation  on  a  Lx L  grid.  During  an  avalanche,  given  that  e  — >  0,  the 
stress  is  not  increased  and  the  relaxation  term  (I)  does  not  participate  in  the  evolution.  This 
approximation  was  performed  for  clarity  purposes,  with  qualitatively  similar  results  with  the 
case  e  =  1.  In  that  case,  the  effect  of  diffusion  is  more  visible  and  avalanches  dissipate  (for 
large  D )  much  smaller  stress  percentage,  since  the  relaxation  term  dominates  the  behavior. 
During  the  avalanche  process,  we  evolve  the  system  by  using  cellular  automata  rules:  when  the 
total  local  stress  crosses  its  aj  threshold,  the  associated  local  slip  0  increases  randomly  with  a 
normal  distribution  with  mean  1  and  variance  1.  In  the  stress  of  Eq.  1,  we  have  also  added  a 
term  for  regularizing  purposes  that  slightly  smoothens  the  slip  profiles.  It  takes  the  form  ctV20 
with  very  small  a  =  0.05.  We  have  checked  for  several  system  sizes  (up  to  642)  that  this 
term  does  not  affect  our  reported  results  in  any  visible  manner.  Also,  we  note  that  this  term 
is  physically  motivated,  in  as  much  as  it  is  connected  to  the  coarse-grained  form  of  the  stress 
generated  by  dislocation  pile-ups  (23).  In  our  simulations  we  used  a  flat  distribution  ranged  in 
(0, 1]  for  the  quenched  disorder  oy(T),  following  a  typical  protocol  (1).  The  kernel  K{ r)  has  a 
continuum  Fourier  representation  K(k )  =  —Ck2k2/(k2  +  k2)2(23),  where  we  set  C  —  1  for 
clarity  purposes  in  our  analysis.  A  modification  of  C  modifies  the  strength  of  disorder  required 
in  the  model  in  order  to  observe  avalanche  behavior,  with  no  other  changes.  In  all  simulations 
using  Eq.  1  (unless  explicitly  mentioned  otherwise)  the  hardening  coefficient  k  is  selected  from 
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the  formula  k  =  2 L°-S5/S0  where  we  chose  S0  =  1000  with  S0  being  approximately  equal  to 
the  cutoff  size  of  the  distribution  that  is  derived  for  D  =  0.  The  reason  for  this  choice  has  to  do 
with  the  fact  that  the  nature  of  the  kernel  (23)  is  such  that  a  fixed  local  hardening  coefficient  k 
does  not  set  the  cutoff  for  the  size  distribution.  Rather,  it  allows  for  a  weak  increase  with  the 
system  size.  However,  for  our  purpose  (studies  of  D  >  0)  it  was  crucial  to  have  well  controlled 
critical  distributions  for  D  —  0,  independent  of  the  system  size,  to  identify  the  concrete  effects 
of  the  relaxation  on  the  distributions.  In  all  plots  we  refer  to  the  value  of  R  =  D/c. 

In  the  simulations  of  Eq.  4,  the  stochastic  equation  was  solved  using  random  variables  that 
follow  the  required  power-law  distribution  with  exponential  cutoff,  generated  with  the  standard 
rejection  method.  During  solving  Eq.  4  numerically,  p  can  jump  above  1,  a  regime  that  we  do 
not  consider.  There  are  several  options  to  deal  with  the  boundary  condition  at  p  =  1  which  are 
numerically  very  similar  for  large  S  and  small  q.  After  a  jump  which  takes  p  >  1, :  i)  p  is  reset 
to  a  random  value  between  0  and  1  (this  method  was  used  for  the  generation  of  Fig.  3  Center), 
ii)  p  is  reset  to  a  specific  value  (for  example,  0),  iii)  p  is  set  to  be  equal  with  the  long-time 
average  1  —  1/S2.  We  shall  reiterate  that  these  crossings  (p  >  1)  are  finite-size  effects  and  do 
not  define  the  system’s  behavior  at  long  times  and  in  the  limit  of  S'  — >  oo,  as  we  verified  in 
both  simulations  of  Eq.  1  (showing  that  the  distribution  “bump”  consistently  vanish  with  the 
system  size)  and  Eq.  4  (showing  that  different  treatments  of  the  p  =  1  boundary  lead  to  the 
same  conclusion  and  phase  boundary  q  ~  1/S). 

The  codes  used  for  the  simulations  of  Eqs.  1  and  4  are  written  in  Python  and  are  available 
upon  request. 

(B)  SOM  Text 

In  this  supporting  online  material  we 
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Figure  SI  -  Duration  Distributions  and  Average  Size  vs.  Durations.  In  associa¬ 
tion  with  the  size  distributions  of  Fig.  1(b),  the  duration  distributions  and  average  size 
for  given  duration  plots,  show  a  consistent  behavior  with  theory  (the  prediction  being 
(S)  ~  T1/Wz  and  at  mean-field  l/crvz  =  2)  (cf.  Fig.  S8).  However,  it  is  worth  be¬ 
ing  sceptical  about  the  importance  of  durations,  given  that  they  are  very  long  compared 
to  the  expectations,  given  the  magnitude  of  typical  individual  dislocation  velocities  (15, 
19).  Here,  we  just  show  that  the  whole  behavior  is  roughly  consistent  with  theoretical 
expectations,  suggesting  that  the  observed  durations  reflect,  perhaps  in  a  serendipitous 
manner,  the  theoretically  expected  behavior. 


(a)  Show  that  the  experimental  results  on  the  durations  of  the  avalanches  are  consistent  with 
the  results  expected  from  our  theory,  and  present  more  details  about  the  methodology  for 
studying  the  rate  dependent  quasi-periodicity. 

(b)  Show  that  the  simulations  of  Eq.  1  are  consistent  with  the  simple  general  model  of  Eq.  4 
which  explains  the  reported  novel  theoretical  behavior. 

(c)  Explore  analytically  the  properties  of  Eq.  4  and  show  that  the  crossover  between  the  two 
found  qualitative  behaviors  can  be  estimated  quite  accurately  using  the  concept  of  the  sur¬ 
vival  probability, 

(d)  Present  more  details  from  the  simulations  we  performed  using  the  model  of  Eq.  1 ,  discuss 
results  from  other  kernels  and,  present  our  conclusions  on  the  generality  of  the  effect. 

(a)  Consistency  of  durations  and  more  details  on  the  experimental  time- 
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Figure  S2  -  Stress  vs.  Strain  and  Strain  vs.  Time  (a)  Selected  stress  -  strain  curves  at 
different  rates,  displaying  similar  hardening  coefficients  and  yield  stresses.  In  the  inset, 
the  strain  vs.  time  behavior  is  shown,  vividly  illustrating  that  the  strain  increases  on 
average  10  and  100  times  slower  in  the  10~5  and  10-6  tests. 

series 

The  study  of  avalanche  durations  in  microplasticity  has  been  scarce  because  the  timescales 
observed  (due  to  a  very  limited  time  resolution  and  apparatus’  sensitivity)  are  many  orders  of 
magnitude  larger  (milliseconds)  than  the  ones  where  the  actual  events  are  actually  taking  place 
(nanoseconds).  Even  though  the  observed  events  are  believed  to  be  well  separated  in  time, 
avoiding  avalanche  overlaps,  the  internal  avalanche  temporal  structure  is  questionable  (2,  3,  4, 
10,  15,  21).  However,  it  is  possible  that  the  observed  large  events  actually  follow  the  statistics 
of  the  dislocation  dynamics’  critical  behavior.  Our  study  of  the  event  durations  and  their  clear 
relation  to  the  event  sizes  suggest  that  this  could  be  the  case:  In  Figure  SI  (a),  we  see  that  the 
duration  distribution  displays  a  strong  rate  dependence  with  the  associated  critical  exponent 
drifting  from  ~  2  to  ~  3,  while  the  average  size  (S)  scales  with  the  avalanche  duration  in  a 
visibly  similar  manner,  having  an  exponent  r^i  2  .  «S>  r\j  T 2)  This  behavior  is  consistent  with 
our  simulations,  as  discussed  in  (c)  of  this  supplement.  The  critical  exponent  for  durations  a 
drifts  from  2  to  2  +  1  =  3,  while  the  relation  between  the  average  size  to  its  duration  remains 
invariant.  In  this  sense,  the  results  resemble  the  effect  of  a  finite  driving  rate  c  (avalanche 
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Figure  S3  -  Windowed  event  timeseries.  Here,  we  display  the  full  timeseries  of  Fig¬ 
ure  1(c)  for  three  selected  samples  studied,  in  order  to  show  that  the  periodicity  we  iden¬ 
tified  is  not  an  artifact  but  characterizes  the  whole  deformation  process.  The  actual  event 
timeseries  (blue  color,  before  smoothening)  is  also  present  in  order  to  show  that  the  peri¬ 
ods  we  identified  are  strongly  correlated  with  strong  peaks  in  the  event  sizes  at  the  high 
points  of  the  periods. 


overlap)  in  the  mean-field  avalanche  model  (3),  wherein  the  average  size  scales  as  rY1/avz  and 
auz  remains  invariant  when  c  is  modified. 

The  phenomenological  hardening  coefficient  used  in  the  simulations  of  Eq.  1  is  held  fixed 
without  assuming  any  time  or  rate  dependence.  However,  experimentally,  we  observe  a  weak 
but  visibly  complicated  rate  dependence  of  the  hardening  coefficient.  The  stress  vs.  strain 
curves  for  three  of  the  samples  studied  (one  at  each  nominal  strain  rate)  are  shown  in  Figure  S2. 
For  these,  the  hardening  coefficient  corresponds  to  the  average  slope  of  each  of  the  curves.  It  is 
clear  that  as  the  rate  becomes  smaller,  the  hardening  coefficient  presents  an  oscillatory  behavior 
with  no  apparent  monotonic  rate  dependence  (with  about  10%  slope  fluctuations).  These  fluctu¬ 
ations  could  certainly  be  attributed  to  the  sample-to-sample  fluctuations.  We  shall  mention  that 
this  fluctuation  takes  place  at  a  time/strain- scale  very  large  compared  to  the  oscillation  scale  we 
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discussed  in  this  paper,  and  thus  is  not  related  to  the  concrete  results  we  presented  in  the  main 
text  and  Figure  1  and  2.  Moreover,  our  main  conclusions  persist  if  a  slight  time  dependence  is 
incorporated  to  the  parameter  k,  resembling  the  experiments. 

Our  microcrystals  are  large  enough  (~  25  //m)  that  their  plastic  deformation  is  representa¬ 
tive  of  a  comparable  subregion  of  the  bulk  (which  has  size  independent  characteristics).  We 
believe  our  observed  quasiperiodic  oscillations  and  the  competition  between  glide  avalanches 
and  slow  relaxation  are  also  happening  in  bulk  plasticity  at  these  deformation  rates.  We  show 
the  full  timeseries  for  the  three  samples  chosen  to  make  Figure  2(a)  in  Figure  S3,  together  with 
the  actual  event  timeseries  before  applying  our  smoothening  procedure.  Using  the  smoothened 
timeseries  as  a  guide  to  the  eye,  one  may  observe  that  these  oscilations  of  the  smoothened 
timeseries  are  related  to  bursts  of  activity  taking  place  at  quasi-periodic  intervals.  In  compar¬ 
ison  to  the  theoretical  plot  of  Figure  2(b)  one  may  also  notice  that  the  experimental  behavior 
has  additional  complexity.  Beyond  the  universal  quasi-periodicity,  there  is  a  non-trivial  event 
structure  during  the  bursts  of  activity,  with  some  “shock-aftershock”  behavior  seen:  However, 
we  have  not  been  able  to  identify  a  detectable  event  asymmetry  by  explicitly  calculating  the 
average  skewness  of  the  burst  activity  in  a  quasiperiod.  It  is  worth  noting  that  he  asymmetry 
in  the  model  of  Eq.  1  vanishes  as  k  is  increasing.  Intuitively,  this  behavior  can  be  understood 
by  considering  the  variance  of  the  p  fluctuations  in  Eq.  4:  When  p0  is  very  close  to  the  parent 
critical  point  (equivalently  k  — >  0),  an  increasing  sequence  of  “shocks”  clearly  dominates  the 
system  behavior. 

(b)  Mutual  consistency  of  the  models  in  Eqs.  1,  4  and  verification  of  the 
minimal  picture 

One  may  argue  that  the  association  of  the  models  of  Eq.  4  and  Eq.  1  is  less  than  obvious 
and  requires  concrete  justification.  Here,  we  provide  robust  numerical  evidence  that  our  inter- 
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(a)  (b) 

Figure  S4  -  The  effect  of  c  on  p  in  the  simulations  of  Eq.  1.  It  would  be  a  naive 
expectation  that  the  increase  of  the  driving  rate  c  is  similar  to  increasing  R.  For  R  =  0, 
but  for  a  finite  c,  it  is  clear  that  (a)  there  is  no  effect  of  c  in-between  the  avalanche  events 
and,  (b)  the  difference  A p/ p  for  the  change  of  p  just  after  an  event  shows  a  dependence 
on  c,  but  it  is  on  the  opposite  direction:  namely,  increasing  c  leads  to  an  increase  of  p 
after  an  avalanche  before  it  starts  decreasing  for  very  large  events.  This  is  a  result  of 
the  fact  that  finite  c  leads  to  overlapping  avalanches,  creating  positive  correlation  among 
events.  Note  that  the  error  bars  display  the  saturated  sample  variation,  while  the  lines 
are  Bezier  interpolants,  weighed  by  the  error  bars.  A  large  error  bar  is  a  signature  that 
there  is  a  non-trivial  intrinsic  saturated  distribution  of  A p,  as  one  would  expect.  Linear 
binning  is  used. 


pretation  of  the  phenomenon  described  by  Eq.  1  is  represented  in  terms  of  Eq.  4.  Our  results 
are  for  simulations  of  Eq.  1  using  the  parameters  mentioned  in  (A)  but  also  with  a  simplified 
infinite -range  kernel  (A'(r  —  r')  =  const.),  which  is  the  mean-field  version  of  Eq.  1.  The  actual 
kernel  used  in  Eq.  1  shows  similar  qualitatively  results,  but  we  choose  to  not  discuss  them  here 
since  the  numerical  correspondence  between  Eqs.  1  and  4  is  qualitative,  but  not  as  quantitative 
and  transparent  (for  the  system  sizes  used).  We  believe  that  the  correspondence  between  Eqs. 
1  and  4  becomes  exact  only  in  the  infinite-range  case  but  the  qualitative  picture  shall  be  right 
for  any  kernel,  even  the  short-range  ones.  Since  the  kernel  used  in  Eq.  1  relevant  for  plasticity 
is  in  the  mean-field  universality  class  (27),  we  believe  that  the  correspondence  of  Eqs.  1  and  4 
is  exact  in  the  infinite-system  size  and  quasistatic  (c  — >  0)  limits.  Also,  we  focus  on  a  small 
L  —  32  system  (except  otherwise  mentioned),  since  the  effects  we  describe  in  this  section  con- 
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(a)  (b) 

Figure  S5  -  The  effect  of  k  on  p  and  reasoning  behind  the  effective  model,  (a)  A  clear 
linear  dependence  of  Ap/p  on  k  is  observed  in  the  simulations  of  Eq.  1 .  Larger  k  denotes 
a  sharper  decrease  of  A p/ p  just  after  an  event,  suggesting  that  pt+1  =  p,  —  /I  k  St ,  where 
A  is  a  model  dependent  constant,  (b)  Moreover,  in  this  model,  one  easily  observes  that 
the  average  event  size  scales  with  k  in  a  simple  manner  (S)  ~  1/k  (see  also  (24)),  and 
holds  for  all  system  sizes.  In  this  way,  we  conclude  that  pt+1  —  pt  —  A^-. 

verge  for  very  small  system  sizes.  Note  that  we  do  not  have,  in  this  paper,  any  analytical  means 
of  “proving”  this  correspondence,  even  though  we  believe  it  is  a  tractable  but  tedious  possibility. 


In  order  to  compare  the  structure  of  Eq.  4  with  Eq.  1,  we  calculate  numerically  the  suscep¬ 
tibility  to  jumping  at  every  time-step,  without  counting  the  steps  inside  an  avalanche.  This  is 
done  by  creating  a  sorted  list  of  the  system’s  local  stresses  cr,  subtracting  the  quenched  thresh¬ 
olds  a f  and,  by  estimating  the  numerical  derivative  at  the  maximum  end  of  this  list  by  averaging 
it  over  a  small  window.  Figure  S4  shows,  as  a  consistency  check,  that  an  increase  of  the  driv¬ 
ing  rate  c  (for  R  =  0  (i.e.  zero  diffusion)  and  small  k)  does  not  affect  p  in  a  visible  manner. 
Between  avalanches  p  remains  constant  as  expected  and,  therefore,  5p(t)  for  time-t  intervals 
is  zero  (Figure  S4(a)).  Also,  after  an  avalanche,  (again  for  R=0  and  small  k)  the  change  in  p 
is  negligible  for  small  event  sizes  S  but,  becomes  significant  for  events  that  are  close  to  the 
system  size;  however,  it  is  not  consistently  decreasing  in  a  monotonic  way.  Rather,  it  increases 
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(a)  (b)  (c) 

Figure  S6  -  The  effect  of  R  in  the  simulations  and  motivation  behind  the  effective 
modeling,  (a)  In  between  events,  p  increases  monotonically  with  a  slope  that  becomes 
larger  with  R.  This  behavior  suggests  that  pt+1  =  pt  +  cd  where  cd  ~  R.  (b),  (c)  After 
an  event,  p  monotonically  decreases  in  a  way  that  is  proportional  to  R  in  both  small  (c) 
and  large  (b)  event  sizes,  but  with  different  proportionality  coefficients.  In  the  effective 
model  we  degenerate  this  difference  by  considering  pt+l  —  pt  —  cdSt/S,  where  here  S 
serves  as  a  constant.  Overall,  our  observations  of  Fig.  5,  6  and  this  one  leads  to  our 
conclusion  that  a  minimal  model  is  pt+i  =  pt  +  cd{  1  —  S/S)  where  cd  ~  R  and  S  is  the 
fixed  point  value  for  the  average  size  at  long  timescales. 


for  intermediate  sizes  while  it  decreases  for  larger  ones;  a  behavior  which  becomes  more  evi¬ 
dent  as  the  rate  becomes  larger  (for  example  notice  c  =  0.1  in  Figure  S4(b)).  The  effect  of  the 
other  free  parameter  of  Eq.  1,  the  hardening  coefficient  k,  is  more  significant  (keeping  R  =  0): 
For  large  k,  the  change  of  p  due  to  an  avalanche  is  simple  5p(S)  ~  —kS,  with  small  finite 
system  size  corrections  (cf.  Figure  S5(a)).  Also,  notice  that  the  saturated  simple  standard  de¬ 
viation  increases  with  the  size,  indicating  the  existence  of  a  corresponding  scaling  distribution. 
Moreover,  if  we  also  utilize  the  fact  that  (5')~  1/k  {{24)  and  Figure  S5(b)),  we  can  conclude 
that  5 p  after  avalanche  events  depends  crucially  on  the  system’s  distance  from  the  critical  point 
A  p/p  ~  S/S  where  S  ={S)  for  R  =  0. 

The  situation  changes  drastically  when  R  is  non-zero:  In-between  avalanches,  p  increases 
in  an  almost  linear  manner  having  a  slope  that  is  directly  proportional  to  R{ cf.  Figure  S6(a)). 
The  result  of  Figure  S6(a)  justifies  the  use  of  the  constant  rate  cd  in  Eq.  4  and,  we  identified  it 
as  a  coarse-grained  relaxation  {i.e.  effectively  proportional  to  R) .  Moreover,  the  presence  of  a 
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Figure  S7  -  A  posteriori  verification  of  the  minimal  picture,  (a)  The  timeseries  of 
the  event  sizes  is  displayed  for  simulations  of  different  Rs,  showing  the  emerging  large 
events  that  dominate  at  large  R.  (b)  After  averaging  the  timeseries  in  (a)  over  a  small 
running  time  window  (50),  histograms  of  the  observed  average  sizes  display  the  expected 
behavior,  as  being  predicted  by  our  minimal  model  (and  picture)  where  the  behavior  of 
p  can  account  for  the  relevant  effect  of  R  on  the  distributions.  Note  that  average  size  is 
related  to  p  with  a  robust  relationship,  even  though  (S) is  much  more  intuitive.  In  the 
picture,  the  sharp  distribution  for  R  =  0  becomes  flatter  as  R  increases. 


finite  R  does  not  significantly  alter  the  response  of  p  to  avalanche  events  for  fixed  k  (cf.  Figure 
S6(b),  (c)),  and  it  ultimately  provides  a  confirmation  of  the  use  of  Eq.  4  as  our  minimal  model. 
We  note  that  the  construction  of  Eq.  4  includes  an  assertion  that  we  neglect  all  other  effects 
the  relaxation  might  have  on  the  internal  structure  of  the  avalanches.  It  is  clear  that  this  is  an 
approximation  that  might  be  worrisome  in  cases  where  the  stress  relaxation  takes  special  forms 
and,  if  it  has  symmetries  that  are  not  included  in  the  interaction  kernel  that  is  activated  during 
avalanches. 

Finally,  we  shall  mention  that  our  construction  is  also  verified  by  the  study  of  the  histograms 
of  p  and  ( S )  (cf.  Figure  S7)  as  R  is  increased.  Even  though  the  structure  of  the  whole  timeseries 
is  altered  in  a  rather  complicated  manner  as  R  increases  (Figure  S7(a)),  the  expected  behavior 
is  clearly  present:  the  observed  ( S )  ~  p  e  (where  9  =  1/2  at  mean-field)  is  spread  over  a  large 
range,  losing  its  well  peaked  structure. 
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(a)  (b)  (c) 

Figure  S8  -  Observables  in  the  simulations  and  consistency,  (a)  The  duration  distri¬ 
butions  are  shown  for  the  simulations  used  for  Fig.  1(a)  and  (b),  showing  that  durations 
display  a  strong  rate  dependence,  as  in  the  experiments  (Fig.  4).  (b)  The  average  size 
vs.  duration  plot,  however,  does  not  display  an  observable  dependence  on  the  rate  (also 
as  in  the  experiments  (Fig.  4)),  apart  from  the  emerging  large  event  flat  profile  at  very 
large  sizes.  The  flat  profile  is  a  signature  of  the  renewing  events  that  saturate  the  cutoff 
behavior  of  the  size  distributions  (bump)  which  are  absent  in  the  thermodynamic  limit. 
This  is  just  a  finite-size  effect.  In  terms  of  a  theoretical  interpretation,  the  absence  of  such 
a  dependence  is  supportive  of  the  idea  that  the  effect  does  not  change  the  crvz  exponent 
of  the  depinning  universality  class,  a  vital  signature  of  distribution  integration,  (c)  In 
the  simulations,  the  increase  of  the  rate  c  (at  fixed  R)  or  the  increase  of  R  (at  fixed  c), 
show  very  similar  behavior.  The  most  characteristic  feature  of  it  though,  is  a  stick-slip 
emergence  in  the  displacement  vs.  time  behavior  that  is  displayed.  As  R  increases,  there 
is  an  emerging  stick-slip  and  quasi-periodic  behavior. 


(c)  Analytical  properties  of  Eq.  4  and  derivation  of  the  crossover  line  q  ~ 

a/S. 

The  probability  distribution  for  Eq.  4  has  the  peculiar  feature  that  it  is  divergent  on  short 
scales,  i.e.  it  is  ultraviolet  divergent.  If  we  take 

P(S)  =  l/AfS-3/2e~s/So  (1) 

where  So  —  1/(1  —  p)  from  the  main  text  and  J\f  a  normalization  constant,  then 

/»oo 

M  =  /  P(S)dS  =  2 /Va  -  2^ti/ S0  +  2v/a/S'0  +  O(l/S03/2)  (2) 
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(a)  (b) 


Figure  S9  -  Probability  distributions  as  the  hardening  coefficient  increases.  Clearly, 
the  exponent  of  the  critical  behavior  remains  unmodified  as  the  hardening  coefficient  in¬ 
creases  and  the  cutoff  moves  to  smaller  sizes.  The  data  are  collapsed  with  S  =  Sk _L5, 
apart  from  the  cutoff  bump  which  fades  away,  consistent  with  the  experiment  and  our 
intuitive  effective  model  (Eq.  3),  where  the  bump  is  not  visible.  The  effective  model 
agrees  with  this  picture,  since  increasing  k  signifies  the  effective  decrease  of  S  and  there¬ 
fore  the  large  increase  of  the  fluctuations  in  p.  However,  in  this  case,  the  approach  to 
p  =  1  becomes  improbable. 


Then,  the  average  avalanche  size  comes  out  to  be, 

poo 

(S)  =  /  SP(S)dS  =  (av^F/2) \fSoj~a  +  a(v r/2  -  1)  +  a/ 2(tt3/2  -  3 vV^o  +  O(l/S0)  (3) 


and  in  the  limit  where  S0  — >  oo  (close  to  the  critical  point)  we  keep  the  leading  term, 


(S)  = 


(4) 


The  fact  that  (S)  does  not  scale  with  S0  but  with  its  square  root,  is  in  the  core  of  the  peculiar 
properties  of  such  special  power-law  distributions:  Generally,  let  our  power  law  be  r  (it  is  3/2 
in  the  case  we  considered),  and  let  us  approximate  an  exponential  cutoff  at  S0  as  a  sharp  cutoff. 
Then  ( S )  is  a  ratio  of  two  integrals  of  power  laws,  and  gives 


(2  ~  t)(Sq~t  —  a2~T) 
(l-rXS^-a1-) 


(5) 
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If  r  <  1 ,  then  we  can  take  a  — >  0  and  S  r^/  So-  If  r  >  2,  then  we  can  take  Sq  — »  oo  and  S  ~  a. 
But  if  1  <  r  <  2,  we  can  only  take  a  — >  0  in  the  numerator  and  S0  — >  oo  in  the  denominator, 
and  get 

S~(t-  l)r(2  -  r)aT_1S'g_T  (6) 


which  agrees  with  our  previous  answer  if  r  =  3/2. 


From  Eq.  4,  the  long-time  behavior  satisfies  the  fixed  point  solution  S  =  (S)  giving  that  p 
will  fluctuate  around  the  value  that  satisfies, 


—(S/a)2 

71 


(7) 


Since  we  are  focused  on  the  calculation  of  the  crossover  line,  it  is  clear  that  it  will  take  place 
when  large  negative  jumps  of  p  are  proliferating.  That  requires  S  being  large  enough  that 
S/S  ~  1  / cfi  leading  to  A pn  ~  —1.  As  mentioned  in  the  main  text,  this  is  roughly  equivalent  to 
saying  that  So/ S  ~  1/c,/.  In  more  detail, 


-Plarge  (  Sm  ) 


P(S)dS 

n 

2e-Sm/S0  2x/7rErfc(ffm/ \/ SpSm) 


S0  \J  Sm/Of 


\J  So/a 


which  in  the  limit  where  Sm 

fjarge 


2 a  l  -  vSr^/S^/a  +  So/a 

cx)  (equivalent  to  saying  that  1  / q,  S0 


(Sp/a)2 

( Sm/a )3/2 


,-Sm/So 


1  —  *^o/ a  +  So/a 

Then,  for  having  A pn  =  —1,  the  number  of  steps  to  wait  is, 


oo), 


(8) 


(9) 


T  = 

1  W 


r^j 


l/Plarge(S/cd) 

Ti  {k  +  4(S'/a)2 


—  2n  (S'/a)) 

(^/a)4 


7T  1 

e 4  cds'/a 


1  7T  £cds/a 

~d2^Sj 


a 
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After  each  such  big  avalanche,  p  spends  a  long  time  ascending  towards  the  critical  point,  which 
corresponds  to  an  integration  procedure  which  alters  the  critical  exponents.  The  number  of  time 
steps  of  integration  (assuming  the  large  jump  gives  p  =  0  just  after  the  jump)  is, 

Ta=~  (11) 

Cd 

neglecting  all  corrections  that  involve  S/S.  The  crossover  line  takes  place  when  Tw/Ta  ~  1, 
which  ultimately  gives  us  the  solution,  in  the  limit  S  — >  oo  that  the  crossover  line  takes  place 
when, 


Cd  = 


3.751 

S/a 


(12) 


The  simulations  which  led  to  Fig.  3  of  the  main  text  are  consistent  with  q  =  8 /(S/a)  as  one 
can  easily  check.  The  discrepancy  of  the  constant  is  due  to  deficiencies  of  the  theoretical  calcu¬ 
lation,  since  it  does  not  include  jumps  smaller  than  1  which  are  abundant  near  the  crossover  line 
and  it  does  not  take  into  account  the  boundary  condition  at  p  =  1  (in  the  simulations,  p  is  reset  to 
a  random  position  after  such  a  jump  contributing  also  to  the  integration).  However,  the  trend  is 
captured  correctly  (~  1  /{S/a)).  The  fact  that  Cd  is  dimensionless  is  consistent  with  all  our  pre¬ 
vious  findings  that  point  to  that  we  are  describing  a  true,  novel  critical  behavior  (quasiperiodic 
oscillating  behavior)  which  is  system-size  independent  (but  certainly,  UV-dependent).  Since 
it  is  a  true,  novel  critical  behavior,  it  is  rather  possible  that  the  two  regimes  are  separated  by 
another  phase  transition  (instead  of  a  crossover),  a  possibility  that  requires  further  investigation 
since  it  is  rather  complicated  to  find  evidence  for. 


(d)  Details  and  conclusions  of  the  simulations  of  Eq.  1 

Here,  we  provide  a  more  complete  picture  of  the  simulations  of  Eq.  1  that  we  described  in 
the  main  text.  While  there  are  several  aspects  that  will  be  carefully  studied  in  subsequent  pub¬ 
lications,  there  are  some  main  observations  and  conclusions  we  can  draw  herein.  In  Figure  S8, 
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we  show  the  behavior  of  durations  for  the  same  simulations  discussed  in  the  main  text  (Figure 
2).  First,  we  note  that  durations  correspond  to  the  times  taken  for  jumps  in  the  displacement 
(cf.  Figure  S8(a)).  As  we  described  in  the  caption  of  Figure  2,  in  order  to  compare  with  ex¬ 
periments  we  needed  to  assume  that  this  duration  is  finite  and  corresponds  to  the  number  of 
cellular  automata  steps  inside  the  avalanche.  This  picture  is  a  posteriori  consistent,  since  the 
experimentally  relevant  time  step  is  significantly  smaller  than  the  timestep  of  the  external  stress 
increase.  However,  this  is  just  a  working  assumption  about  the  dynamics  of  the  avalanches  and 
the  situation  could  certainly  be  more  complex.  Overall,  the  durations  display  distributions  that 
are  consistent  with  the  experiments  and  also  consistent  with  our  minimal  picture.  As  seen  in 
Figure  S8(b),  the  durations  have  a  scaling  distribution  and  drifting  exponent  a  from  2.0  ±0.1 
to  3.0  ±  0.2,  with  changing  strain  rate.  This  finding  is  consistent  with  Figure  SI.  The  behavior 
of  (S)  also  displays  a  power-law  with  respect  to  the  avalanche  duration  (( S )  r\ j  Tl/(yvz)  which 
stays  visibly  invariant  as  c  decreases  (cf.  Figure  S8(c)).  There  is  a  plateau  for  large  durations, 
we  have  checked  that  this  feature  disappears  as  the  system  size  increases. 

In  Figure  S9,  we  show  the  behavior  of  the  size  distributions  (for  fixed  relaxation  and  driving 
rates  R  =  2,  c  =  0.000125)  as  k  increases  for  L  =  128.  For  k  =  0.01  the  system  already 
displays  strong  quasi-periodic  behavior  that  we  described.  A  fixed  system  size  controls  the 
avalanche  size  cutoff  scale  ,S0  according  to  the  formula  S0  ~  2L°-85/ k.  From  Figure  S9,  it  is 
straightforward  to  observe  that  the  critical  exponents  of  sizes  (Figure  S9(a))  and  durations  (Fig¬ 
ure  S9(b))  are  invariant  to  the  change  in  the  cutoff  scale  and,  therefore,  we  conclude  that  the 
distance  from  the  parent  critical  point  (controlled  by  k )  does  not  alter  the  integration  effect  (see 
Eq.  4)  on  the  critical  distributions.  This  result,  together  with  the  inset  of  Figure  3(a),  signifies 
that  the  integration  effect  we  report  does  not  depend  on  the  distance  of  the  system  from  the 
critical  point  before  R  is  turned  on.  However,  it  is  also  easy  to  notice  that  the  distance  from 
the  critical  point  strongly  affects  the  functional  form  of  the  distribution  near  the  cutoff.  When 
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the  system  for  R  —  0  is  near  to  the  parent  critical  point  (small  k)  the  size  and  duration  distri¬ 
butions  display  a  characteristic  “bump”  near  the  cutoff.  This  bump  is  the  finite-size  signature 
of  avalanches  that  percolate,  which  equivalently  corresponds  to  p  becoming  equal  to  1  before 
these  avalanches  take  place.  However,  as  expected  this  bump  vanishes  when  k  increases,  but 
with  no  visible  change  in  the  critical  exponents.  In  a  qualitatively  similar  manner,  the  bump 
vanishes  when  the  system  size  increases  (Figure  S  10(b))  or  when  R  decreases  (Figure  S  10(a)). 
Overall,  this  evidence  signifies  that  the  phenomenon  we  describe  corresponds  to  true  critical 
behavior  as  relaxation  processes  are  present  in  the  system,  and  not  a  special  finite  size  effect. 

Finally,  we  present  some  relevant  conclusions  about  the  simulations  that  we  shall  describe  in 
detail  in  a  future  publication:  First,  the  quasi-periodicity  and  associated  integration  (see  Eq.  4 
and  discussion  therein)  is  present  for  all  interaction  kernels  we  checked  in  two  dimensions: 
K( k)  =  |k|.  k2,  or  const,  all  display  the  reported  phenomenon  as  R  increases.  The  kernel  |k| 
was  studied  in  (20)  for  a  single  slip  plane  of  mixed  dislocations,  and  k2  corresponds  to  a  local 
kernel  which  leads  to  avalanche  critical  behavior  in  the  short-range  interface  depinning  uni¬ 
versality  class  (7).  Second,  the  quasi-periodicity  and  associated  integration  appears  for  several 
functional  forms  for  the  stress  relaxation  term  of  Eq.  1;  We  have  explicitly  checked  a)  local 
0  diffusion  (~  V20),  b)  non-local  diffusion  (~  ((0)  —  0)),  c)  local  stress  diffusion  (acting 
between  0  slips,  0  evolves  in  a  way  that  da/dt  ~  V2a  is  satisfied).  We  believe  that  the  case 
studied  in  (9)  strongly  resembles  the  much  simpler  model  of  K( k)  =  |k|  and  the  (c)  type  of 
relaxation.  All  of  these  models  belong  to  the  same  class  displaying  the  general  phenomenon 
described  in  this  paper. 
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Figure  S10  -  Probability  distributions  as  R  increases  and  stability  with  system  size. 

(a)  As  R  increases,  the  model  of  Eq.  1  displays  quasiperiodic  oscillations  and  the  expo¬ 
nent  r  of  the  size  distribution  drifts  to  2.5  from  the  original  1.5.  Notice  that  the  distribu¬ 
tion  displays  a  “bump”  at  the  cutoff,  which  grows  as  R  increases.  This  bump  corresponds 
to  large  events  close  to  the  critical  point,  which  cause  large  changes  of  the  jumping  sus¬ 
ceptibility.  However,  the  bump  is  a  finite-size  effect,  since  it  goes  away  when  the  system 
size  increases  or  when  k  decreases,  (b)  For  fixed  R,  the  system  size  dependence  is  stud¬ 
ied  for  P(S),  showing  a  consistent  scaling  collapse  with  S  ~  SL~ 0  94  (inset).  We  shall 
note  that  for  R  —  0,  the  size  distribution  does  not  scale  with  the  system  size  (the  hard¬ 
ening  coefficient  has  been  chosen  accordingly  to  scale  with  the  system  size  in  order  to 
achieve  this),  in  order  to  isolate  the  system  size  dependence  related  to  R.  The  only  non¬ 
collapsing  part  is  the  “bump”  region,  where  the  bump  consistently  disappears  in  the  large 
system-size  limit,  signifying  that  the  system  never  crosses  the  parent  critical  point  to  the 
infinite-avalanche  phase  (i.e.  large  events  are  scale  invariant,  part  of  the  scaling  behavior 
and  do  not  scale  with  the  system  size).  The  details  of  this  behavior  will  be  discussed 
elsewhere. 
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